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the hypothesis of randomness corresponding to a Poisson point process. We 
first compute the exact first and second moment of the Ripley K-statistic 
under the homogeneous Poisson point process model. Then we prove the 
asymptotic normality of a vector of such statistics for different scales and 
compute its covariance matrix. From these results, we derive a test statistic 
that is chi-square distributed. By a Monte-Carlo study, we check that the 
test is numerically tractable even for large data sets and also correct when 
only a hundred of points are observed. 
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1. Introduction 



Analysis of point patterns is relevant in many sciences: cell biology, ecology or 
spatial economics. The observation of clusters in point locations is considered as 
a hint for non observable dynamics. For example the clustering of tree locations 
in a forest may come from better soil conditions or from spreading of seeds of a 
same mature individual; but clusters are also observed in random distribution 
as a Poisson point process sample. It is therefore essential to distinguish be- 
tween clusters resulting from relevant interactions or from complete randomness. 
Ripley (1976, 1977) — is a widely used tool to quantify the structure of point 
patterns, especially in ecology, and is well referenced in handbooks (Ripley, 
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1981; Diggle, 1983; Stoyan et al., 1987; Cressie, 1993; M0ller & Waagepetersen, 
2004; Ilian et al., 2008). Up to a renormalization by the intensity of the process, 
this statistic denoted here K(r) estimates the expectation K(r) of the number 
of neighbors at distance less than r of a point in the sample. The observed K(r) 
is compared to the value of K(r) for a homogeneous Poisson point process with 
the same intensity as the data, chosen as a null hypothesis: the Poisson point 
process is characterized by an independence of point locations, modelling an 
absence of interactions between individuals in ecosystems. In this case K(r) is 
simply the mean number of points in a ball of radius r divided by the intensity, 
that is Trr 2 . If K(r) is significantly larger than irr 2 (respectively smaller), the 
process is considered as aggregated (respectively over-dispersed) at distance r. 
To decide if the difference is statistically significant, we build a test of the 
Poisson process hypothesis; we need to know the distribution of K{r) for this 
process. But even the variance is not known and statistical methods generally 
rely on Monte-Carlo simulations. Ripley (1979) used them to get confidence 
intervals. Starting from previous results (Saunders & Funk, 1977), he also gave 
critical values for the L function, a normalized version of K introduced by 
Besag (1977). These critical values are valid asymptotically, for a large number 
of points but low intensity, so that both edge effects and point-pair dependence 
can be neglected. Further computations of confidence interval bands based on 
simulation have been proposed in Koen (1991) and corrected in Chiu (2007). 
But the simulation is a practical issue for large point patterns, because com- 
putation time is roughly proportional to the square of the number of points 
(one has to calculate the distances between all pairs of points) multiplied by the 
number of simulations. 

We propose here to compute the exact variance of the Ripley statistic. Ward & Ferrandino 
(1999) studied this variance. But they ignored that point pairs are not inde- 
pendent even though points are (eq. A8, p. 235), thus their derivation of the 
variance of K{r) was erroneous. The right way to compute the covariance is to 
consider that it is a [/-statistic as remarked in Ripley (1979), then to use the 
Hoffding decomposition. As the variance is not enough to build a test, we study 
the distribution of the statistic. We prove its asymptotic normality as the size 
of the observation window grows. It is then easy to build an asymptotically 
Gaussian test. 

Another concern is to test simultaneously the aggregation/dispersion at dif- 
ferent scales. This is rarely correctly achieved in practical computations with 
Monte-Carlo simulations. The confidence bands or test rejection zone are often 
determined without taking the dependence between the numbers of neighbors 
at different scales into account. As an exception Duranton & Overman (2005) 
provide a heuristic multiscale test. In our main theorem, we consider a set of 
scales (ri, . . . , r^), compute the covariance matrix of the K(ri) and prove the 
asymptotic normality for the vector (K(rx), . . . , K{r<i)). From this we propose 
the first rigorous multiscale test of randomness for point patterns. 
The paper is built as follows: Section 2 introduces the precise definition of 
K{r) and the current definition of K{r). In Section 3, after the definition of 
our statistics (no edge-effects correction, known or unknown intensity), we list 
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the main results of the paper: exact bias due to the edge effects and exact 
variance of K(r) for a homogeneous Poisson process with known or unknown 
intensity; covariance between K{r) and K(r') for two different distances r and 
r' . The main theorem contains the convergence of the vector (K(r\), . . . , K(rd)) 
to a Gaussian distribution with explicit covariance in the following asymptotic 
framework: data from the same process are collected on growing squares of ob- 
servation. These results allow a simple, multiscale and efficient test procedure of 
the Poisson process hypothesis. Section 4 provides a Monte-Carlo study of the 
test and Section 5 gives our conclusions. The last section contains the proofs. 
Technical integration lemmas are postponed in the appendix. 

2. Definition of the Ripley Jf-function 

We recall the characterizations of the dependence of the locations for a general 
point process X over R 2 . We refer to the presentation of M0ller & Waagepetersen 
(2004). 

2.1. Definitions 

For a point process X, define the point process X^ on R 2 x R 2 of all the cou- 
ples of two different points of the original process. The intensity of this new 
process gives information on the simultaneous presence of points in the original 
process. Denote p( 2 \x,y) its density (called the second-order product density). 
The Poisson process of density p(x) is such that p^ 2 \x,y) = p(x)p(y). 
The Ripley statistic is a way to estimate the density p^ 2 \x,y). Precisely it is 
an estimate of the integral on test sets of the ratio g(x,y) = p^ (x,y)/ p(x)p(y). 
The function g(x, y) characterizes the fact that the points x and y appear simul- 
taneously in the samples of X. If g(x, y) = 1, the points appear independently. 
If g(x,y) < 1, they tend to exclude each other; if g(x,y) > 1, they appear more 
frequently together. 

We assume the translation invariance of the point process: g(x,y) = g(x — y). 
In order to estimate the function g, we define its integral as the set function K.. 
Let A be a Borel set: 

K(A) = [ g(x)dx. 

J A 

If we also assume that the point process is isotropic, we define the Ripley K- 
function as 

K(r) =fC(B(x,r)), 

where B(x,r) is the closed ball with center x and radius r. The translation 
invariance implies that K.(B(x,r)) does not depend on x. For example, if the 
process is a Poisson process then g{x) = 1 and K{r) = irr 2 . We define the 
Rpley statistic that estimates the fiT-function. Let A be a bounded Borel set of 
the plane R 2 , m the Lebesgue measure and p an estimator of the local intensity 
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of the process; for a realization S of the point process X, S = {X\, . . . , Xjv}, 
the Ripley statistic is defined by 



3. Main results 



This section presents the theoretical results on the Ripley statistic and the 
resulting test. 



3. 1 . Definitions 

Throughout the paper, we refer to the indicator function I, the expectation 
e rj „, the centred indicator function h and its conditional expectation h\. We 
gather here these definitions. 

Let n be an integer; A n denotes the square [0, n] 2 ; U is a random location 
in A n with an uniform random distribution; its density is 1/n 2 with respect 
to the Lebesgue measure d^id^ over A n . V is a random location with the 
same distribution as U and independent of U. We denote d(x, y) the Euclidean 
distance between x and y in the plane, and I{A} the indicator function of set 
A. We define e r>n = E( l{d(U,V) < r}), h(x,y,r) = I{d(x, y) < r} - e r> „ and 
hx{x,r) =E(h(U,V,r)\ V = x). 

3.2. Assumptions 

We assume that X is a homogeneous Poisson process on R 2 with intensity p. 
We consider that the data are available on the square A n . S = {Xi, . . . , X^} 
is the sample of observed points. We consider two cases: 

1. If the intensity p is known, the Ripley statistic is expressed as 

M0 = i E I{d(**,*j)<r}. 

2. If the intensity p is unknown, we choose to estimate p 2 by the unbiased 
estimator p 2 = N(N - l)/n 4 (Stoyan & Stoyan, 2000) and define 

^.»W = M^iiy E I{«,Jfi)<r}. 



3.3. Bias 



It is known that a large number of neighbors of the points located near the 
edges of A n may lie outside A n causing a bias in the estimation. We compute 
the bias due to this edge effect. 
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Proposition 1. Assume that r/n < 1/2. 

EK U r)-K(r) = r 2 (_| + ^) . 
~ / 8r r 2 

-'''""'(•- 1 + ^)! ,+ <"''"'"') ■ 

Notes: 

• The assumption that r/n is less than 1/2 means that at least some balls 
of radius r are included in the square A n . 

• The additional term for K^^n corresponds to the probability to draw a 
sample with zero or one point in the square. This probability is so low 
that the term gives a zero contribution as soon as the mean number of 
points pn 2 is larger than 20. 

• The proof may be adapted for a convex polygon of perimeter Ln to com- 
pute the first order term of the bias; for u = 1 or 2: 

^ or„2 

EK u , n (r)-K(r) = + 

6 n 




3-4- Variance 

We compute the covariance matrix of K UiTl (r) for it = 1 or 2. We get an exact 
computation for the variance, that can be used for any value of n. 

Proposition 2. For < r < r' , 



war(£i, n (r)) = ^ + + ^-Efc?(tf,r), 



P* P P 

cov(K ltn (r),K lin (r')) = gg^ + ^'"^ + ^-cov^{JJ^),h x {U,r)), 



P 2 P P 



var(K 2>n (r)) = 2n 4 



( I{N > 1} 
\N(N - 1) 



Am ( I{N>l}(N-2) \ Bu2nr 
+ 4nE { N(N-l) ) Eh ^ 

+ n 4 e- pn2 (1 + pn 2 ) (l - e""" 2 - pn 2 e~P n ^ e 2 >n , 
co-y(A' 2! „(r),^ 2 .n(r')) = 2n 4 E f - ^ > ^ (e r ,» - e r /,„e r , n ) 



\N(N- 1) 

+ 4n4g ( /{ ^-l')~ 2) ) «»"(fai(^0,fti(^r)) 
+ n^""™ 2 (1 + pn 2 ) fl- e-"" 2 - pn 2 e-P n2 ^j 
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where 

irr 2 8r 3 r 4 

er - n ~ l^~~3n3 + 2rt [, 

^ l2/TT . r 5 /8 256 \ r 6 /ll 56 \ 8 r 7 lr 8 

E/if(t/,r) = — -7T + ~f — 7T M =■ -. 

1V ' ' n 5 \3 45 / n 6 \48 9 / 3 n 7 4 n 8 

./Voies: 

• The variances of both estimators are exact and can be computed at any 
precision, as inverse moments of the Poisson variable correspond to fast 
converging series. But these series may be difficult to evaluate with math- 
ematical softwares, because of the large value of the Poisson parameter. 

• The covariances are not explicit because the terms cov (hf(U, r'), hf (U, r)) 
involve terms that have to be numerically integrated. 

• The leading terms of the variances of K\ >n (r) and K2, n (r) as n tends to 
infinity are 2nr 2 /n 2 p 2 + Anr A /n 2 p and 2nr 2 /n 2 p 2 . 



3.5. Central Limit Theorem 



We show that a normalized vector of Ripley statistics for different r converges in 
distribution to a normal vector. Let W(0, S) denote the Gaussian multivariate 
centred distribution with covariance matrix E. 

Theorem 1. Let d be an integer, < r\ < . . . < r<j a set of reals and define 
IC u ,n = (K u ,n(ri), • • • , K Utn (rd)). Then n^/p(JC u . n - ir(r\, . . . r 2 d )) converges in 
distribution to 7V(0, E) as n tends to infinity, where for s and t in {1, . . . , d} 

• ifu=l, E Sit = h 4tt r s r t . 

P 

2n(r 2 A r 2 ) 
P 

Note: The first term of the variance corresponds to a situation where the 
couples of points are independent from each others; this was used as an ap- 
proximation without proof in Ward & Ferrandino (1999); our work proves that 
the actual variance and limit process are different in the first case and that the 
approximation holds only in the second case. 

3.6. Applications to test statistics 

From Theorem 1, we deduce that T u = E _1 / 2 /C Mi „ is asymptotically Af(0,Id) 
distributed. For the hypothesis 

Ho: X is a homogeneous Poisson process of intensity p 

we use T 2 = ||T U ||2 as a test statistic with rejection zone for the level a: 

T 2 > X 2 a (d). 
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Figure 1. Comparison of normalized variances for Ki(l), p = 5 



where Xa(d) is the (1 — a)-quantile of the x 2 (d) distribution. 
Note: the covariance matrix £ depends on the intensity parameter p, so that 
in the case of the unknown parameter we have to use an estimate of p in the 
formula defining E. 

4. Simulations 

We study the empirical variance of the proposed statistics by a Monte-Carlo 
simulation. Then we apply the test procedure to simulated data sets, observe 
the number of rejections and compare it to the level of the test. 

4-1- Variance 

We simulate a sample of 1000 repetitions with p = 5 and compare (after renor- 
malization by n^fp) the empirical variance and the exact computed variance 
with the limit variance for different value of n (figure 1). With 1000 repetitions, 
the oscillations of the empirical variance are still large; we will use a larger num- 
ber of repetitions in the following study of the test. 

The convergence of the computed variance to the limit value is not so fast and 
for applications with hundreds of points (corresponding in figure 1 to n < 15) 
the distance between the variances is still large. A preliminary study, not pre- 
sented here, showed that the test procedure is perturbed by an small error in the 
covariance matrix, as we tried simplified versions of the covariance by bounding 
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Table 1 

Percentile of rejection over 10000 repetitions of the test with level a = 0.05. 



Poisson 


T* 


T> 


Ti 


T* 


T 2 


n 


= 30 


p = 


1 


r = 


(1,2,5) 


5.40 


5.04 


5.20 


5.01 


5.10 


n 


= 10 


p = 


5 


r = 


(1,2,5) 


5.61* 


5.40 


5.19 


5.38 


5.37 


n 


= 10 


p = 


5 


r = 


(1,2, ...,10) 


5.13 


5.32 


5.76* 


6.67* 


6.01* 


n 


= 10 


p = 


1 


r = 


(1,2,5) 


5.67* 


5.86* 


5.81* 


5.30 


5.25 


n 


= 10 


p = 


.5 


r = 


(1,2,5) 


5.52* 


5.73* 


5.52* 


5.60* 


4.91 


n 


= 10 


p = 


.2 


r = 


(1,2,5) 


6.40* 


6.84* 


6.59* 


6.59* 


5.22 



or ignoring the corner contribution C'(A^ 3 ) (see in the proof section). It is cru- 
cial to use an accurate computation of the covariance matrix to have a correct 
approximation of the square root inverse matrix £ -1 / 2 . Therefore we will use 
the exact formula instead of the asymptotic formula in the test procedure. 



4.2. Test 

In the known parameter case, the computation of the test statistic T\ is straight- 
forward; we also build a statistic T* using the empirical covariance matrix of 
the sample. The advantage of T* is that it is orthogonal by construction and 
should lead to better results. But the covariance matrix is not observable when 
we dispose of one sample, so that the test procedure based on T* is unfeasible. It 
is an idealized version, used to compare the corresponding number of rejections. 
To avoid the statistical dependence between the sample and the estimator of 
the covariance matrix, we also build a statistic T[ where we generate a addi- 
tional independent sample of the Poisson process with intensity p to compute 
the empirical covariance matrix. 

In the unknown parameter case, the computation of the test statistic T2 is simi- 
lar. In the variance formula the unknown parameter p is replaced by the estima- 
tor N/n 2 . We also choose to replace the expectation E ( I{iV > 1}/(N(N - 1))) 
by the observed value 1/(N(N-1)) and E ( 1{N > 1}(N - 2)/(N(N - 1))) by 
(N — 2)/(N(N — 1)), because the dispersion of a Poisson variable is low with 
respect to the expectation when its intensity is large. The construction of T 2 * is 
the same as for T*. The case of T'^ is not studied because, as p is unknown, one 
would have to generate an additional sample for each estimated value of p. 
The test output is a Bernoulli random variable with parameter a. With a suf- 
ficient index of repetition m, the mean number of rejection is close to a normal 
variable with expectation a and variance a(l — a)/ m. We consider that the test 
works when the observed frequency of rejection is in the 95% Gaussian confi- 
dence interval [a — 1.96-y/a(l — a)/m, a + 1.96 x /a(l - a)/m]. With m = 10000 
and a = 0.05, the interval is [0.0457; 0.0543] so that the percentile of rejection in 
table 1 should lie in [4.57; 5.43]. Stars indicate the values outside the confidence 
interval. 

The performances in the case of a known parameter (Ti, T* and T[) are good 
except when the number of points is small. The unfeasible tests T* and T[ based 
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Table 2 

Percentile of rejection over 10000 repetitions of the test with level a = 0.05. 



Thomas 


T 2 


n = 10 («, ft, a) = (1, 5, 3) r = (l,2,5) 
n = 10 (k, h, a) = (0.5, 10, 0.5) r = (1, 2, 5) 


71.6 
100 



on the empirical covariance have no better performance than the test T\. The 
error of the empirical covariance is probably still to large. The only exception is 
the third line where a large number of values of r are considered simultaneously. 
The test T2 performs better than T\ for small data sets. The only exception is 
the case of a large number of scales. The poor performance of T\ and T2 in this 
case may result from numerical instabilities in the covariance matrix inversion 
as its dimension is larger. The departure from normality may also be larger in 
this case (some classes of inter-point distances being weakly represented in the 
sample). With this exception, the test based on T2 works perfectly. 
In table 2, we investigate the power of the test T2 by simulating two Thomas 
cluster processes (Thomas, 1949). A Thomas process is a Neyman-Scott process; 
the germs of the clusters are drawn as a sample of a homogeneous Poisson process 
of intensity n . For each germ, an inhomogeneous Poisson process is drawn with 
intensity measure /i/, where / is the density of the Gaussian two-dimensional 
vector centered on the germ and with independent coordinates of variance a. 
The Thomas process results from the superposition of these Poisson processes. 
The germs are not conserved. The parameters of the two processes are such that 
clusters are not visually detectable in the first process and evident in the second 
one. The test rejects 71% of the first sample and systematically the second 
one. The test is more powerful than a visual observation of the data, detecting 
invisible clusters. A rigorous analysis of the distribution of the statistic for 
dependent point process models should allow to conclude on the power of our 
test but such a study is beyond the scope of this paper. 

5. Conclusion 

We provide an efficient test of the null hypothesis of a homogeneous Poisson 
process for point patterns in a square domain. This is a theoretical and practical 
improvement on preexisting methods: Monte-Carlo simulations are untractable 
when the number of points increases. With a personal computer, calculating K 
for 10,000 simulations of a 10,000-point set is not feasible (or it will take months). 
Marcon & Puech (2003) applied K to a, 36,000-point data set (the largest ever 
published as far as we know), but had to limit the number of simulations to 20. 
We suggest to change the treatment of edge effects. Instead of correcting edge 
effect on each sample to reduce the bias, we compute the exact bias. The use 
of sample correction (for each point of the data) has not been questioned since 
Ripley's original paper, except by Ward & Ferrandino (1999). 
We also point out that the test can be used on samples with a few dozens of 
points as encountered in actual data sets. It works correctly with such small 



G. Lang and E. Marcon/ Testing randomness of point patterns 



10 



data sets, even if it is based on asymptotic normality. This is due to the fact 
that the bias and variance are known exactly and not asymptotically; the non- 
normality of the statistics for small data sets seems to have lesser effects than 
approximating the variance. 

Our work should be extended in two directions: to other domain shapes that 
are of interest for the practitioners and to 3-dimensional data for high resolution 
medical imagery. A further study of the asymptotics of the distribution of K (r) 
for dependent point process models such as Markov or Cox processes should 
also be achieved to inform on the power of our test. 

6. Proofs 

6. 1 . Proof of proposition 1 

Recall that U and V are two independent uniform variables on A n . The expec- 
tations of the Ripley statistics are 



E 




E(N(N-1)) 



E( l{d(U,V) < r}) 




n 2 P(7V>l)E( \{d(U, V) < r}) 




The following lemma allows to conclude: 





(1) 
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Figure 2. Zones in the square 



where (see figure 2) 

• (interior) ={£, £ is at distance larger than r from the boundary} 

• (edge) A\ ={£, £ is at distance less than r from an edge, larger than r 
from the others} 

• (two edges) A\ ={£, £ is at distance less than r from two edges and larger 
than r from the corner} 

• (corner) A^ ={£, £ is at distance less than r from the corner} 

Note that A^, A^ and A\ are composed of four parts that contribute identi- 
cally. We establish formulas only for one of these parts. 



Lemma 2. Define function g(x) = arccos(a:) — x\l — x 2 . 

/ I{d{£,i]) < r}dn = nr 2 . 
JneAn 

If£€ A 2 n , with n - r < £i < n, xi = £(n - £i), 

I{d(£,n) < r}dn = r 2 (ir - g{x{)) 



s 



'r/GAr, 

If £ & A^, with n — r < £i < n, n — r < £2 < n and (xi,x 2 ) = -(n — £1, n — £2), 
I{d(£,ri) < r}dn = r 2 (ji - g(xf) - g(x 2 )). 



J 



If£E A\, with n - r < £1 < n, n - r < £ 2 < n and (xi,x 2 ) = £(n - £i,n — 

g(x x )+g(x 2 y 



[ I{d{£,n) <r}d V = r 2 [^+ Xl x 2 

JrieA n \ 4 
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Figure 3. Geometrical interpretation of g 



Note: Function g(x) is the area of the part of a ball of radius 1 that lies outside 
the square when the ball intersects one of its edges (see figure 3). 
Proof. For the interior points £ £ A\, B(£, r) C A n . 
Let £ £ A\. We compute the area of r) n A n . 

/ rf) < r}drj = + 2r 2 / y/l - t 2 dt 

JneA n 2 J Q 



it — arccos(a;i ) + x\ y 1 — : 
= r A (tt - g( Xl )) . 

Note that r 2 g{x) is the part of the ball that lies out of the square A n if the 
center is at distance xr from the edge of the square. 

Let £ £ A^. Here the ball intersects two edges of the square and the area of 
B(£,r)nA n is 

rj) < r}di] = r 2 (tt - g(iEi) - fif(x 2 )) . 

Let ^ G Af t . Divide the ball into four quarters along axes parallel to the coor- 
dinate axes. One of the quarter is inside the square, two intersect the edges, 
leaving outside an area equal to (g{xi) + g(x2))/2. The area of the intersection 
of the last quarter with the square is x\X2 so that the area of r) PI A n is 

J a L{d{£,r))<r}dr) = r \^—+x 1 x 2 J. □ 

Proof of lemma 1 (continued) . The left-hand side of (1) is m(A} l )TTr 2 = tt(u — 
2r) 2 r 2 . Recall that A 2 n is composed of four parts that contribute identically. 
We integrate function g. 
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Lemma 3. 



G{x) = I g(u)du = xarccos(x) - y/l-x 2 + - (1- x 2 ) 3/2 + -. 



Proof. Changing variables and integrating by parts 

/>arccos(x) 

arccosfw)^ = — / tsm(t)dt 

Jir/2 



/>arccos(:c) 

[tcos(i)]™ (a;) + / cos(t)dt 



tt/2 



x arccos 



(x) - y/l - x 2 + 1. 



Changing the variable w = vl — u 2 , we get 



y uy/l-u 2 du = J v 2 dv = 3 ((! - ^ 2 ) 3/2 ~ 1 



□ 



Then the contribution (2) is equal to 

^ 2 i! r2< < K ~ 9 ^ dx = 4r 3 (n-2r)(vr-G(l)) = ^4tt - r 3 (n-2r). 

We consider A 3 ; the domain of integration is symmetric in [x\,X2) so that the 
contribution (3) is equal to 



4r 4 ^ dxi ^ (■n-2g{x l ))dx 2 =<ir i ^tt (l - ^) - 2 ^ g(xi)dxij ^ dx 2 



From Lemma 6, 

5(^1)^1 / dx 2 =G(l)- / g(xi) Jl - x\dx x = - - — . 

J^/T^l Jo v 3 lb 

We consider ^4^; the contribution (4) is equal to 



so that contribution (3) is equal to r [ 47T 

2 3 



4r 4 y dxi J y— + £ix 2 - £(xi)Jcte 2 = r 4 ( — + - - 4 y g(xi)y 1 - x^cixi 

= r 4 

Gathering the four contributions, we get 

r '™ n 2 ^ ( n ) ( 3) n ( n ) ( 6 ) n 2 ^ 

r 2 / _8r lr^\ 
n 2 \ 3 n 2 n 2 
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6. 2. Proof of proposition 2 

We decompose the variance of K Sj a„ i r ) by conditioning the variable with respect 
to the number N of points in the sample. Conditionally to N, K S: A n {r) has 
the form of a [/-statistic. Then we apply the Hoffding decomposition to this 
[/-statistic. 

For s = 1, 2, we use the relation 

var(K s , An (r)) = v&vE(K SiAn (r)\N) +Evax(K a , An (r)\N). 
We first consider the conditional expectation of K Si A n {r). 

E{K x , n {r)\N) = JL I J2 E l{d(X t ,X 3 ) <r})= N ( N ~ 1 ) e r,n ^ 
n z p z \ / n z p z 



,2 



N 



E(K 2<n (r)\N) = " Yl E WW*' u i) ^ r > = " 2 ^« ^ > X >" 

Because iV is a Poisson variable with intensity pn 2 

EN 2 (N-1) 2 = EN(N - l)(N — 2)(N - 3) 

+ 4EN(N-l)(N-2) + 2EN(N-l) 
= p i n 8 +4p 3 n 6 + 2p 2 n\ 

var N{N~- 1) = 4/rV + 2p 2 n 4 . (5) 

Then 

wt& / Mart (4pn 2 + 2)e 2 
varE(Jf lin (r)|iVj = ! -. (6) 

varE(i? 2 ,n(r)|iV) = n 4 P{N > 1}(1 - ¥{N > l})e 2 „ 

= n 4 e-"" 2 (1 + pn 2 ) (l - e^" 2 (l + pn 2 )) e 2 „. (7) 

We compute the conditional variances. 

var(^ lin (r)|iV) = ^var ( £ h^X^r) ) , 
71 P \i&=i J 

4 ( N 

vav(K 2 , n (r)\N) = —-—— var £ h{X. ll X j 

Conditionally to iV, the locations of the points are independent and uniformly 
distributed variables [7, over A n . We introduce the Hoffding decomposition of 
the [/-statistic kernel h: 

h{x,y,r) = hi(x,r) + fti(j/,r) + h 2 (x,y,r), 
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where h\[x) = E(h(U,V,r)\V = x), (U,V) being two independent uniform 
random variables on A„ . 

Then E1n(U,r) = and E(h 2 (U,V,r)\U) = E(h 2 {U,V,r)\V) = 0, so that 

var h(U, V, r) = var hi(U, r) + var hi(V, r) + var h 2 iJJ , V, r) 

= ZEhl(U,r)+vzrh 2 (U,V,r). 

From 

N N N 

h(U i ,U j ,r)=2{N-l)J2hi(U i ,r)+ £ h 2 (U i ,U j ,r). 

i^j=l i=l i#J = l 



we get 

vax{K x>n {r)\N) = 



4(JV- 1) 



-^-var (j2 h i(Ui,r)] + ^var ( £ ha^U^r)] 



~4 4 

no n p 



n p n p 
4jV(jV-l)(jV-2) ^ 2 ^ 2JV(JV-1) 
= nV Eftl(y,r) + nV Var/l(C/ ' F ' r) ' 

Now var h(U, V, r) = e r .„ — e 2 n and using factorial moments of the Poisson 
distribution 

~ 4n 2 2 

E var(X 1 .„(r)|iV) = Ehf(U,r) + — (e r ,„ - e 2 ) . (8) 

p p 2 

Lemma 4 gives the exact value of E/i 2 ([/, r). With relations (6) and (8), we get 

(«■ ( \\ 2e r , n 4n 2 e 2 ,, t 4^^,^^ 

var(Ai n (r)) = — 5 — I 1 EftJ ({/,-, r) 

P 



P 2 


1 


p 




'27rr 2 


47r 2 r 4 


I? ( 


v P 2 






'16r 3 


/32tt 


( 




+ (- 






/59tt 


^ ( 




U + 



45 
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Similarly 



1{N > 1}{N - 2) 2 ^ : 2n 4 I{7V>1} 



var(JT 2 , n (r)|iV) = jv(JV -I) + iV(iV - 1) varfe (^ F ' r )' 

Evar(/? 2 ,„(r)|iV) = E ( ^ ^(t^.r) 



From this and relation (7), we get 

var(/? 2 ,„(r)) = 2n 4 E ~ <n) 

+ nV 5 " 2 (l + pn 2 ) (l - e-"" 2 - pnV" 2 ) e 2 



■T. n * 



We now apply the same decomposition to cov (i£i )Tl (r), iTi )n (r')), 

cov(E(5- 1 , n (r')|JV),E(X- 1 , n (r)|i\0) = ^!±^£^ , ( 9 ) 



cov (£!,„(/), ^i,„(r)|iV) = 4(iV 4 4 1)2 cov (Vfti^yj.V/ii^r) 



/ JV JV \ 

—cov ^ h2(Ui,Uj,r'), J2 MUi,U jt r)\ 



rrp 
2N(N-1) 



n 4 p 4 



cov (h(U,V,r'),h(U,V,r)). 



- - 4n 2 2 

E cov(#i,„(r'), lfi in (r)|jV) = cov (fti(J7, r'), fei (C/, r)) + ^(e r! „ - e r > n e r n ). 

P P 

To compute cov (h\(U, r'), h\(U, r)), the square A n should now be split into 16 
different zones according to the 4 zones of the preceding section with respect to 
r and the 4 zones with respect to r'. Because of inclusions, the actual number 
of zones to consider is reduced to 9. The corresponding computation is easy in 
the center zone, but can not be achieved in a close form in the edge bands and 
in the corner. We consider the following zones: 

• (interior) A^ 1 ={£, £ is at distance larger than r' from the boundary}, 
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• (interior-edge) A^ 2 ={£, £ is at distance between r and r' from an edge, 
larger than r' from the others}, 

• (edge) ^4 2 ' 2 ={£, £ is at distance less than r from an edge, larger than r' 
from the others}, 

• (corner) A 3 ^ 3 ={£, £ is at distance less than r' from two edges}. 

Denoting x\ — Mn-fi) and x[ = h{n — £1) we get 

hi^.rOM-X^r) = (^" _er '' n )(5" _er ' n ) on 
7rr' 2 ?'' 2 . , A / 7rr 2 



^r' 2 r' 2 



-2fl , ( a; i) — 5" - e r,n 25(^1) on A 



n 2 \ 8r 
7r ^e r 



vr J \ n n 
~ 2n 2n 2 ' 



covihiXjyihiX^r)) = C(A^) + C(A^ 2 ) + C(A 2 n 2 ) + C(A 3 n 3 ) 



2r' x 2 



C(A J n 

f 1r'\ r' 3 r 2 f 1 

C{A\i 2 ) = 4 1 - — —b r ,n / (b r ,, n - gix'^dx', 

\ Tl J n Jr/r' 

C(A 2 /) = 4(l-—) 7 ^-[(b rl>n -g(rx 1 /r'))(b r!n -g(x 1 ))dx 1 . 
\ n / n Jo 

The first integral may be expressed in terms of function G, the second integral 
is elliptic and has to be numerically evaluated; as the integrand is bounded and 
very smooth this can be achieved without difficulties. To compute the term 
C(A 3 1 ' 3 ), we rewrite the different values of function hi with the help of indicator 
functions: 

h A \{x,r) = b r>n l{xi > l;x 2 > 1} 

h A 2{x 1 r) = (b r>n - g(x 2 )) l{xi > 1; x 2 < 1} + (b r>n - g(xi)) \{x 2 > 1; x\ < 1} 
h A3 (x, r) = (b r<n - g{xi) - g{x 2 )) l{x\ < 1; x 2 < 1; x\+x\> 1} 
h Ai (x,r) = (b r<n - 7r/4 + xix 2 - (g(xt) + g(x 2 ))/2) \{x\ + x\ < 1} 

For x' = p-(n-£i,n-£ 2 ) 

2 14, /.l fl 4 4 
71 JO JO i=1 i=3 

and this integral also can be numerically evaluated. 

Note: the whole computation of this term of the covariance could be numerically 
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achieved, but it is preferable to use an exact computation whenever it is possible. 
The case of the covariance of K 2iTl {r) is analogous: 

cov(E(K2, n (r')\N),E(K 2 , n (r)\N)) = n 4 e -"" 2 (l +pn 2 ){l-^-P n \l+pn 2 ))e r , n e r , n 



E cov{k 2 , n {r'),k 2 , n (r)\N) = 4n 4 E ( ^n^-I) ~ ) C0V (M^y),M^)) 



2n 4 E ^ Z~jy ) ( er <™ ~ e r',ne r .n) ■ □ 



6.5. Proof of Theorem 1. 



We show that any linear combination of the Ki^ n (rt) is asymptotically normal. 
Let A = (Ai, ... \d) be a vector of real coefficients. Define Z\ ~ 2t=i ^t^i,n( r t)- 
We use the Bernstein blocks technique (Bernstein, 1939): we divide the square 
A n into squares of side p with p = o(n). These squares are separated by gaps of 
width 2rd so that the sums over couples of points in each square are indepen- 
dent. The couples of points with at least one point in the gaps give a negligible 
contribution, so that the statistic Z\ is equivalent to a sum of independent vari- 
ables and asymptotically normal. 

Set p = n 1 / 4 . Assume that the Euclidean division of n by (p + 2rd) gives 
a quotient a and a remainder q. For / = 0, . . . , a, we define the segment 
h = [(p + 2r d )/, (p + 2r d )l + p — 1]. We order the set {0, . . . , a} 2 by the lexico- 
graphic order. To any integer i such that 1 < i < k = (a + l) 2 , corresponds an 
element (ji, j 2 ) of this set; we define the block P^ n = Ij, xlj 2 and Q = A n \UiPi !n 
the set of points that are in none of the P^n's. For each block and Q, we 
define the partial sums: 

A t l{d{X h X m )<r t }, 

d 

^A t l{d{X u X m )<r t ] 
=i 

t I{d(X h X m )<r t }. 



k 

n^fp{Zi - EZi) = y^(Mi, n - Eu tyn ) + y^X v i,n - E«i,n) + W n - Ew n , 
i=l i=l 

We show that the sum of the Ui_ n converges in distribution to a Gaussian variable 
and that the other term are negligible in 2 . We check the conditions of the 
following CLT adapted from Bardet et al. (2008). 



1 d 

1 Xi^x m eP t , n t=i 



then 



rip 



3/2 



np 



£ 5 

d 



X t ^X m £Q t=l 
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Theorem 2. Let (zi. n )o<i<fc(n) & e an array of random variables satisfying 

1. There exists 5 > such that E\zi^ n \ 2+S tends to as n tends to 
infinity, 

2. Ylj—ri varzi in tends to a 2 as n tends to infinity, 

then y^_n %i,n tends in distribution to Af(0,a 2 ) as n tends to infinity. 

To check Condition 1, we compute the fourth order moment of Ui^ n — Euj i?l . 
Let Ni be the number of points of S that fall in Pi, n . Define 

d d 

f( x , v)=Y Xt ( l{d(x, y) < r t } ~ e r , p ) = Y X t h(x, y, r t ) 



E{{ Ui , n - Eu itn ) 4 \Ni) = -l^E ( Y f( U l' U ™) 
Denote f\ and f 2 the decomposing functions of /: 

E(/i(l/,)) = 0, E(/i(Uj)/ 2 (l/i, Cm)) = E{h{U m )h{U u U m )) = 0, for [7 ; and U„ 
two independent uniform variables on Pi, n - 

Ni Ni Ni 

Y f(Ui,U m ) = 2(N i -l)J2f 1 (U l )+ Y h{U h U m ). 

l=£m=l 1=1 ljtm=l 



Note that |ft,i(x,r)| < nr 2 p 2 so that f\ is bounded by Cp 2 . 

Define M x = E fe^ /i(C7|)) 4 . Then Mi - NiE(J}(U))+6Ni(Ni-i)E(ft(U)) s 



and 



E(iVi - l) 4 Aii = 0(1). 



Define M2 = E \Y^i^m=\ fa{Ui, U m )J . Because f 2 is zero mean with respect 
to one coordinate, only the products where variables appear at least two times 



contribute. 

Mo 



Ni Ni 

8 ®fl(Ui,U m ) + 16 ®fl(Ui,U u )fl(U m ,U u ) 

Ni 

32 Y E fi(Ui,U m )f 2 (U m ,U u )f 2 (U l ,U u ) 

Ni 

4 Y, ^fl(Ui,U m )f 2 (U u ,U v ) 

Ni 

16 Y ^h(UuU m )h{U mi U u )f 2 {U u ,U v )h{U v ,Ui). 
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Because f 2 is bounded, EM 2 = 0(ENi(Ni - 1)(JV, - 2)(N t - 3)) = 0(p 8 ), so 
that 

fc 

i=0 

As p = 71 1 / 4 , we get condition 1. 

To check condition 2, note that the vector (Ki.p i (n), . . . , -ftTi^p^rd)) has a co- 
variance matrix E p defined by Proposition 2 by substituting p to n in the ex- 

pressions. The Uj ; „ = p ^2 t=1 ^t(Ki,p,{ r t) — Elf^ (r t )) are i.i.d variables 

4 

with variance equal to ^-^A'EpA. But p 2 pE p tends to S as p tends to infinity 
and 

Vvaru, „ = -^-A t T.„A — > A*EA 

i=0 

so that X)f=i tends in distribution to A/ r (0,A*SA). 

Note that the are fc independent variables. Denote Ni iri the number of 
points Xi in the boundary region Pi. Td of Pj >n such that the ball B(Xi,rd) 
intersects Q and let D(Xi) denote this intersection. Note that 

EN i>rd = pm(P i>rd ) < C P r d . 



va,rv it „ 

where 



ENi, rd EN%T> 2 {X m e D(Xi)\X m G 0} 
p 3 m(P i)r J(m 2 (Q) + m(Q)) (^|y) = 0(p)- 

E^ I{A m G B(A,,2r d )} I{X„ G n D(X m )} 

I— I m— 1 tt=l 

EiV 2 rd P{X m G S(X;,r d )|X m G P i>rd }EN Q F{X u G £>(*0I*« G Q} 

+ (^y) "«> (5^) - o { p). 



V ;=i TO =1 



/=1 m— 1 it— 1 

N i>rd N t , rd N Q 

E E EE ^ G 

/— 1 m— 1 it— 1 



< 

T 2 = 
< 
< 
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and var (^2i=i v i,nj — O (fcp/n 2 ) = O (p^ 1 ), so that this sum is negligible in 
2 . Similarly 

C ( Nq V C 

var( W „) < — E MT l{X m eB(X h r d )}\ <—(T 1 +T 2 ). 



, Z^m=l 



where 



Ti = E^^ UX m eB(X h r d )} 



1=1 m=l 



< E7V Q (iV Q - l)P{X m G BpQ,r d )|X ro eQ}< m 2 (Q)- 



m(Q) 

T2 = E EEE I{X m GS(X i ,r d )}I{X u GB(X,,r d )} 



/=1 m— 1 u—1 



< EiV 2 (7V Q - l)P 2 {X m G r d )|X m G Q} 

< (m 3 (Q) + 2m 2 (Q)) ' 



Then var (w„) = O (m(Q)/n 2 ) = O (p- 1 ) and w n is negligible in . 
Consider now K2, n(f)- Define Z2 = Y^t=i ^t^2,n(i't) = 4zv,n^i where Ajv, ra = 
wfey- We have E ( A n]u) = 1 and from (5), var (A N ] n ) = — + -j-j. 
For S > 0, the Markov inequality gives 

F(|^-l|>^)<^^- 

Then, with 5 = n" 1 / 4 

00 00 4 2 

e p (i^:« - !i > ™~ 1/4 ) < e ^ + < 00. 

n— 1 n— 1 ' ^ 

From the Borel-Cantelli lemma, we get that Aj J 1 n converges a.s. to 1. By the 
Slutsky lemma, A^, n Z\ converges in distribution to AT(0, A'SA). □ 

6-4- Computation of Eh^(U,r) 
Lemma 4. 

, 2/ . r 5 /8 256\ r 6 /ll 56 \ 8 r 7 1 r 8 
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Proof: From the computation of the bias, denoting x% = Mn — we get 

2 

hi(£,r) = — ^ — e r>n on A\ 



r 



= -j(7r-ff(ii))-e,., n on4 



r 2 



(tt - g{xx) - 5(2:2)) - e r . n on A 3 



-j I — + X!X 2 ) - e r ,„ on A* 

2r\ 2 r 4 , / 2r\ r 5 r 6 , 

Effti Xj.r ))' = tt 2 1-- -j-<„ + 4 1-- _r 1 + 4-(r 2 + r 3 

n ) n \ n J n° n° 

Ti = / (tt - gix^fdx! (10) 
Jo 

T-2 = I dxi / (tt - - g(x 2 )) 2 dx2 (11) 



r, = /W^(f +Ms _£M±lfe2))fe. (12) 

To compute these three terms, we need integral computations on function g. 
Lemma 5. For n> 1, 



I n = / u 2n 1 arccos(u)c?w 



?r(2n) 



n2 2 "+ 2 (n!) 2 ' 

J n = [ ?i 2 Vl - u 2 dw = -(2n + 2)J„+i + 2nl„. 
Jo 



TT 2 


1 


= 16 H 


" 4" 


TT 3 


7T 


= 48^ 


" 4' 


3tt/64, Ji 


= 7T 



(13) 
(14) 



Lemma 6. 



J\(u)y/l-v?du = ^. (15) 



(i 



.9 2 N^ = y-g- (16) 



5 2 ( U )vr^& = -. (l?) 

[ g(u) G (y/l - u 2 ) du = —- — + -. (18) 
Jo V / 96 48 9 v ' 
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Proofs are postponed in the appendix. Using these lemmas, we get 

T X = ^-2^0(1)+ j\ 2 (x 1 )dx 1 ^n 2 -^-^. (19) 

T 2 = 7r 2 ( 1 — - J - An / g{xi)dxx / dx 2 + 2 / g 2 (xi)dxi / dx 2 

v 4/ J hfT^ Jo " Jx/i^ 



+2 / g{xi)dxx [ 

Jo J\/i z xl 



_g{x 2 )dx 2 . 



/■l /■ 1 g 3 

From the computation of the bias, —An / g(x\)dxi I dx-> = 1 . 

Jo " J^/T^ 3 4 



2 / g 2 (x 1 )dx 1 / Gte 2 = 2 / g 2 (xi)cfei-2 / Jl - xjg 2 (x 1 )dx 1 = — — — — . 

Jo J^/T^ Jo Jo v 3 45 24 



From (16), (17) and (18), we get 



2^ g{x 1 )dx 1 J ^_g{x 2 )dx 2 = 2G 2 (l)-2^ 3 (.ti)G ^1 - x 2 ^) da* 



V 

Adding these results, we obtain 



7T 3 57T 

'48 + 24' 



7T 3 9 9tt 128 , x 

T2 = 16 + ~8 45"' ^ 



To compute T3, we write 



— — + x x dxi / x 2 da;2-— / g{xi)dxi / cto 2 
4 Jo Jo 1 Jo Jo 



- / 3 2 (xi)dxi / dx 2 +— I x x dx x I x 2 dx 2 
*> Jo Jo z Jo Jo 

- Ig{xi)dxi I g(x 2 )dx 2 -2 1 xig(xi)da;i / x 2 dx 2 . 
* Jo Jo Jo Jo 



1 rV 1 -*! 1 f 1 I 1 7T 

x\dx x J xldx 2 = -J x\(l - x\)^ 1 - x\dx x = -{J x - J 2 ) = — . 

From (15), — - / g(x{)dxi / dx 2 = — — . 

From (17), - j g 2 (x 1 )dx 1 J dx 2 = - — . 

— J xxdxi J x 2 dx 2 = — J xx(l - x 2 )dx x = —. 



192 ~ 96 + 9' 



1 Z" 1 /V 1 -^ 
From (18), - / g(x 1 )dx 1 / g(x 2 )dx 2 
* Jo Jo 



G. Lang and E. Marcon/ Testing randomness of point patterns 



24 



,1 ry/l-*l /■! 

—2 / a;i(?(a;i)da;i / x^dx-z = / (xf — xi) <?(a:i)<ixi 

Jo Jo Jo 

Adding these results, we get 

7T 3 19tt 2 

T 3 = 1 h -. 

6 16 192 9 

Gathering (19), (20) and (21) gives the result. □ 
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Appendix A: Integration lemmas 
A . 1 . Proof of Lemma 5 

Integrating by parts 



,1 ,tt/2 l «r/2 

/ u 2n ~ 1 aiccos(u)du= I t cos 2n_1 (i) sm{t)dt = — / cos 2n (t)dt. 
Jo Jo 2n Jo 



Using De Moivre formula 

cos 2n (t) = (2cos(2ni) + 2 ( 2 ™ 1 cos(2(n - l)t) 



2n 
n 



Only the last term gives a non zero integral, giving the result for /„ 



(u 2n+2 - u 2n )(~(l - u 2 y 1/2 )du 

I' 1 

l ) arccos(u)] Q / ((n + 2)u 2n+1 - nu 271 ^ 1 ) arccos(u)dM 



,2n+2 2n 



and the term under brackets is zero, giving the result. 



r 1 r^/ 2 r 

/ yl — u 2 aiccos(u)du = / tsin 2 (t)dt = 
Jo Jo Jo 



<"* /2 t tcos(2t) 



dt 



o 
16 
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A. 2. Proof of lemma 6 

Equation (15) follows from equation (13). 

Write g 2 (u) = arccos 2 (u) + u 2 — u 4 — 2uy/l — u 2 arccos(u) and 

arccos 2 (u)du = / t 2 sm(t)dt = - [t 2 cos(t)] +2 tcos(t)dt 
Jo Jo 

7f/2 



2[tsin(t)]o /2 + 2 / sin{t)dt = n - 2, 



/■ 2 , 112 

(u — u )du = = — . 

v ; 3 5 15 



uyl — u 2 arccos(ii)dw = / t cos(t) sin 2 (t)e?i 

t/2 ^ f 7r/2 



zr/2 



n 



^sin 3 (t) 



3 7o 



sin 6 (t)dt 



7T i / ,7r/2 i r /2 

= / sm(t)dt + - / cos 2 (t) sin(t)dt 

6 3 Jo 3 7 

7T 1 1 r 3/,oW2 7T 2 

= fcos 3 (t) L = . 

6 3 9 L WJo 6 9 

Collecting the three parts yields to (16). 



g 2 (u)\/l — u 2 du = j yl — u 2 arccos 2 (u)du 
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-21-- — 
48 ' 8 V 8 64/16 32 

Write G (^J 1 — x 2 ^j = \J \ — x 2 (J^ — arccos(a;)^ + ^ — x + ^ 

g{x)G — £ 2 ) ^ X = J ^ ~ x2 _ arccos ( :r )) arccos(.T)dx 
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— / (x — .t 3 ) ^ — arccos(x)^ dx 
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